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We have performed the first global QCD analysis to include the CCFR and NuTeV 
dimuon data, which provide direct constraints on the strange and anti-strange parton distri- 
butions, s(x) and s(x). To explore the strangeness sector, we adopt a general parametriza- 
tion of the non-perturbative s(x), s(x) functions satisfying basic QCD requirements. We 
find that the strangeness asymmetry, as represented by the momentum integral [S~] = 
J x[s(x) — s(x)]dx, is sensitive to the dimuon data provided the theoretical QCD constraints 
are enforced. We use the Lagrange Multiplier method to probe the quality of the global fit 
as a function of [S~] and find —0.001 < [S~] < 0.004. Representative parton distribution 
sets spanning this range are given. Comparisons with previous work are made. 
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1 Introduction 



The recent measurements of both neutrino and antineutrino production of dimuon final 
states (charm signal) by the CCFR and NuTeV collaborations [1] provide the first promising 
direct experimental constraints on the strange and anti-strange quark distributions of the 
nucleon, s(x) and s(x). In addition to the intrinsic interest in nucleon structure [2], the 
strange asymmetry (s — s) has important implications on the precision measurement of the 
Weinberg angle in deep inelastic scattering of neutrinos. [3-8] We report here the first global 
QCD analysis that includes the new dimuon data, using the methods developed by the CTEQ 
collaboration, specifically to explore the strange and anti-strange parton parameter space. a 
In previous global analyses, information on s and s has resided only in inclusive cross 
sections for neutral and charged current DIS. The reliability of the extraction of the quite 
small s and s components' 3 (from differences of large cross sections measured in different 
experiments) was always in considerable doubt. For this reason, most global fits adopted 
the assumption s(x) = s(x) = k(u + d)/2; this approximation was inferred from the earlier 
combined neutrino and antineutrino dimuon experiments [11] which extract a consistent 
value of 

^_ Jdx x[s(x,Q 2 ) + s(x,Q 2 )] _ Q4 
f dx x [u(x, Q 2 ) + d(x, Q 2 )] 

at some low value of Q. The recent high-statistics dimuon measurements of [1] provide 
greater accuracy, and the high purity of separate neutrino and anti-neutrino events offers the 
first opportunity to study the difference s(x) — s(x). Neutrino induced dimuon production, 
(y/u)N — > fi + fi~X, proceeds primarily through the subprocesses W + s — > c and W~s — > c 
respectively, and hence provides independent information on s and s. 

We present the first global QCD analysis that includes this new dimuon data. The 
new results demonstrate, first of all, that the strangeness asymmetry, as measured by the 
momentum integral 

[S~] = f x[s(x) - s(x)]dx , (2) 
Jo 

is indeed more sensitive to the dimuon data than to the other DIS data. We then use the 
recently developed Lagrange multiplier method of global analysis to explore the range of 
uncertainty of [S~]. In this first report, we concentrate on [S~], the integrated strangeness 
asymmetry, which represents a new parton degree of freedom in the nucleon heretofore largely 
unexplored, and which has immediate impact on precision electroweak physics because of 
the NuTeV anomaly. [3] A full exploration of the density functions s(x) and s(x) will be 
presented later. 

a A preliminary version of this study was reported at the Lepton Photon 2003 International Symposium 
(LP2003), Fermilab, August 2003. Cf. P. Gambino [9] and R. Thome [10], published in the Proceedings of 
LP2003. 

b The strangeness content of the nucleon, as measured by the momentum fraction carried by s or s, is of 
order 3% at Q = 1.5 GcV. 
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We begin by describing the general features of the strangeness sector of the nucleon 
structure in the QCD framework, and our general parametrization of that sector. A brief 
review and discussion of the pQCD calculations that are relevant to the interpretation of 
the dimuon data is then followed by the main results of the global analysis, with emphasis 
on concrete representative global fits relevant for probing the strangeness asymmetry, and 
on Lagrange multiplier results for the integrated momentum asymmetry (2). The paper 
concludes with a summary of the extensive studies performed beyond the examples given, 
comparisons to previous work on strangeness asymmetry, and conclusions. 

2 General properties of s(x) — s(x) and its first two mo- 
ments 

Before we discuss the concrete fitting procedure and the results in the following Sections, 
we feel it is instructive to formulate the qualitative expectations that are based on general 
QCD requirements (before any data). 

2.1 Requirements 

For each Q, let us define the strangeness number densities s ± (x) and their integrals [s ] by 

[s^ = / s ± (x)dx = I [s{x) ± s{x)\ dx , (3) 
Jo Jo 

and the momentum densities 5 ,± (x) and integrals [S^] by 

[S*] = / S ± (x)dx= [ x[s(x)±s(x)]dx . (4) 
Jo Jo 

In the QCD parton model, certain features of these quantities are necessary: 

• The parton distributions s(x,Q 2 ) and s(x,Q 2 ) c [or equivalently s ± (x)}, are parametrized 
at some low (but still perturbative) scale Q ; the full Q-dependence is then determined 
by DGLAP evolution. d 

• The strangeness number sum rule for the nucleon requires 

[s~] = (for all Q). (5) 

A necessary corollary is that the density function s~(x) must be less singular than 1/x as 
x — > for all Q. As DGLAP evolution preserves Eq. (5) (which is a consequence of the 
conservation of the strange vector current = ^ S 7 M ^ S ), it suffices to impose it at Qq. 

c As we had done already in Eqs. (3), (4) above, we will suppress the obvious Q 2 dependence from now 

on. 

d For definiteness, the CTEQ PDFs evolve from Q = m c = 1.3 GeV. 
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• The momentum sum rule requires 

[S + ] = 1 - S (for all Q), (6) 

where S represents the momentum fraction of all non-strange partons. Through this 
condition the global inclusive DIS and other data — which directly affect S — indirectly 
also constrain [S + ]. 

• In the limit x — > (high energy and fixed Q) , Regge considerations and the Pomeranchuk 
theorem predict s~(x)/s + (x) —> 0. 

2.2 Expectations 

From the above general constraints, we draw the following conclusions: 

(i) The number sum rule, Eq. (5), implies that a graph of s~(x) must cross the x-axis at 
least once in the interval < x < 1; and the areas bounded by the curve above and below 
the x-axis must be equal. 

(ii) Assuming a simple scenario in which there is — as supported by theoretical models [2] 
based on Lambda-Kaon fluctuations — only one zero crossing, either s~(x) < in the low 
x region and s~(x) > in the high x region, or vice versa. The two possibilities imply 
[S~] > or [S~] < 0, respectively, because the momentum integral suppresses the small x 
region and enhances the large x region. 

(hi) The low-x behaviour of |s~(x)| ~ x 13 - will then be correlated with the size of l^ - ]!: 
The steeper the function s~(x), the larger ^"(x)! has to be at large x. 

To illustrate point (ii), we preview a representative CTEQ fit (to be discussed in full detail 
in Sec. 4) in Fig. 1 and juxtapose it with results from previous literature. According to the 
parametrizations of s(x) and s(x) used by the CCFR-NuTeV dimuon study [1, 12], s~(x) is 
negative in the x range covered by the experiment (0.01 < x < 0.3). e A previous detailed 
global analysis of inclusive data by Barone et al. [4] (BPZ), finds that s~(x) is positive 
in the large x region. These previous results are shown in Fig. 1 as the dot-dashed and 
dashed curves, respectively. In light of the theoretical constraints discussed above, both 
these results would hint at the first possibility mentioned above, i.e. [S~] > 0. This is not 
the case, however, for the CCFR-NuTeV curve, as it does not cross the x-axis; it violates 
the strangeness number sum rule and cannot be taken to represent the physics outside a 
limited window in x. The BPZ curve does satisfy the sum rule; we note (Fig. 1) it has 
two zero-crossings. The solid curves (Class B) in the two plots of Fig. 1 provide a concrete 
example of s~(x) and S~(x) following the above requirements and expectations. 

Since the preliminary version of this work was reported, the CCFR-NuTeV collaboration has emphasized 
that their most recent analysis favors an integrated strangeness asymmetry that is consistent with zero 
(K. MacFarland and P. Spentzouris, communications at the LP03 Symposium and WIN03 Workshop). More 
definitive studies are needed to clarify the situation. Sec further discussions in Sec. 5. 
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Figure 1: Different functional behaviors of the strangeness number asymmetry function 
s~(x) and momentum asymmetry function S~(x): comparison of our central fit "B" (solid, 
discussed in detail later in the text) with those of BPZ (dashed) and CCFR-NuTeV (dot- 
dashed). The horizontal axis is linear in z = x 1 ^ 3 so that both large and small x regions are 
adequately represented; the functions are multiplied by a Jacobian factor dx/dz so that the 
area under the curve is the corresponding integral over x. 



To illustrate point (iii), we preview various classes of solutions (again, to be discussed 
in detail below) in the upper plot of Fig. 2. Because the experimental constraints are weak 
or non-existent in the very small x region, say x < 0.01, the detailed behavior of s~(x) is 
unconstrained in this region. However, this uncertainty at small x is considerably reduced 
in S~(x), as demonstrated by the curves of the lower plot. Thus the above observations 
concerning [S~] are affected only mildly by the uncertainty of the very small x behavior, 
unless that behavior is so extreme that the small x region provides a significant contribution 
to the number sum rule. We will refrain from exploiting such a mathematical possibility as 
long as it does not seem to be motivated by any physics. f 

f A residual bias of the results on the assumed functional forms, or rather on the rejection of some forms 
as merely mathematical and not physical, is unavoidable in parton analysis. One can imagine the extreme 
(unphysical) scenario of spike- like structures escaping "detection" in the x — ► region but, nevertheless, 
affecting sum rules. 
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Figure 2: Typical strangeness asymmetry s (x) and the associated momentum asymmetry 
S~(x), as obtained in our global analysis. The axes are the same as in Fig. 1. 

3 General parametrization of the strangeness distribu- 
tions 

To explore the strangeness sector of the parton structure of the nucleon, we need a suit- 
able parametrization of s(x) and s(x) (or equivalently s ± (x)) at a fixed scale Qq. This 
parametrization must satisfy the theoretical requirements specified above, and it should be 
as general as possible so that the allowed functional space can be fully explored. A general 
form is essential, so that our conclusions are not artifacts of the parametrization, but truly 
reflect the experimental and theoretical constraints. In the following we explain our choice of 
such a parametrization. Full details including explicit parameters are given in an Appendix. 

It is more natural to parametrize the s ± (x, Qq) functions independently (rather than 
s and s) since they satisfy different QCD evolution equations: pure non-singlet for s~ and 
mixed singlet /non-singlet for s + . We use the following parametrizations, 

s+^Qo) = A x M {l-x) M P + {x-A^A^...) (7) 
s~(x,Qo) = s + (x, Qq) tanh[ax 6 (l — x) c P-(x; xq, d, e, ...)] (8) 

where P+{x; A 3 , . . . ) is a positive definite, smooth function in the interval (0, 1), depending 
on additional parameters A3,... such as are used for u,d,g,... in most CTEQ [13] and 



6 



other global analyses [14, 15]; and 



P_(s) = 




where the crossing point xq is determined by the strangeness number sum rule [s~] = 0, and 
the parameters d, e, . . . are optional, depending on how much detail is accessible with the 
existing constraints. Important features of this parametrization are the following: 

• The strangeness quantum number sum rule, [s~] = 0, is satisfied by the choice of Xq. The 
parameter xq has a physical interpretation: it is the "crossing point" where s~(x) = 0. (If 
d and/or e are not zero, there can be additional zeros of s~(x); in practice — as explained 
in the previous section and below in Sec. 4.2 — we restrict our analysis to solutions with a 
single crossing only.) 

• The fact that the tanh function has absolute value less than 1 ensures positivity of s(x) 
and s(x). The fact that tanh is a monotonic function guarantees that the function s~(x) 
can be made as general as necessary by the choice of P_(x). 

• The small- X behavior of s~(x) must be such that the integral [s~] converges (before the 
root xo is determined). Let j3- = A\ + b; then Eq. (8) implies 

s~(x) ~ x^~ as x — > 0. (10) 

The convergence of [s~] is guaranteed if /3_ > —1, i.e. the parameter b is chosen in the 
range b > —1 — A±. 

Because P±(x) can be made as general as necessary, the choice in Eqs. (7)-(9) is capable of 
exploring the full strangeness parameter space allowed by data in the PQCD framework. 

Detailed formulas for all flavors used in the actual analysis are presented in the Ap- 
pendix. 

4 Global Analysis 

We now describe the global QCD analysis, which includes all relevant experimental data 
and implements the theoretical ideas outlined above. This may be considered an extension 
of the on-going CTEQ program of global analysis. Several new elements (compared to the 
latest CTEQ6M [13] analysis) are present. On the experimental side, we have added the 
CDHSW inclusive F 2 and F 3 data sets [16], and the CCFR-NuTeV dimuon data sets [1]. On 
the theoretical side, we have expanded the parameter space to include the strangeness sector 
as discussed in Sec. 3. 

Compared to the global analyses of BPZ [4], which also allow s, the major difference 
experimentally is our inclusion of the dimuon data, which provide a direct handle on s and 
s; and, theoretically, the generality and naturalness of our parametrization of the strange 
distributions. 6 Since the results of [4] on strangeness asymmetry rely on small differences 
g Ref. [4] parametrizes s(x) and s(x) rather than s ± (a;). 
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of inclusive DIS charged-current and neutral-current measurements, BPZ performed the 
analysis at the cross section level, applying uniform procedures to treat data from different 
experiments in the comparison to theory. Considering small differences between inclusive 
cross sections, the possible strange asymmetry is but one of many sources that could lead to 
such differences. 

As the dimuon data more directly constrain the strange PDFs, this is an important 
new element to our fit. For all inclusive DIS processes we use the standard procedure of 
comparing theory with the published F2 and F3 structure function data. In our analysis, 
the fit to charged-current (neutrino) inclusive structure functions is dominated by the high 
statistics CCFR data. Although we have included the earlier inclusive CDHSW data (which 
play a prominent role in the analysis of [4]), they have no discernible influence on the results 
presented below. 

To include the CCFR-NuTeV neutrino and antineutrino dimuon production data in a 
global QCD analysis is not a straightforward task. The experimental measurement is pre- 
sented as a series of "forward differential cross sections" with kinematic cuts, whereas the 
theoretical quantities that are most directly related to the parton distribution analysis are 
the underlying (semi-) inclusive "charm quark production cross sections". The gap between 
the two is bridged using a Monte Carlo program that incorporates kinematic cuts as well as 
fragmentation and decay models. In our analysis, we use a Pythia program provided by the 
CCFR-NuTeV collaboration to do this efficiency-correction. 11 This Monte Carlo calculation 
is done in the spirit and the framework of leading-order (LO) QCD. CTEQ5L parton distri- 
butions and Peterson fragmentation functions were used. The parameters of the model were 
tuned to reproduce, as closely as possible, the detailed differential dimuon cross sections 
published in [1]. 



4.1 da vN ^ cX in QCD 

At LO in PQCD, the cross section formula for vN — > cX is [17] 

- 2 G 2 fMnEu \ v <*\ dxdy 



■2 



is'^Q 2 ) + 0{a s ) , (11) 



2M N E„C 

with the CKM matrix element |V^ S | and where the Barnett-Gottschalk parameter [18, 19] 



2 

m 



e s '\} + <?) (12) 

h We thank Tim Bolton and Max Goncharov, in particular, for providing this program, as well as assistance 
in its use. Their help was vital for carrying out this project. 
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approaches Bjorken-x as Q — > oo (relative to m c = 1.3 GeV). The quantity s'(£, Q 2 ) e ff in 
Eq. (11) includes Cabibbo suppressed contributions in neutrino scattering via 

W I 2 

s' = s + ^^d (13) 

\v cd \ 2 

with obvious adjustments for the anti-neutrino case 1 . 

The NLO corrections to £s'(£, Q 2 ) c s, defined via the perturbative series 

i 

are generically of the form 

e,'(e,Q 2 )2cx J2 f® H f ( 15 ) 

f=g,s> 

with (g) denoting a convolution integral over parton momentum. These were first calculated 
more than 20 years ago [19]. Later calculations [17, 20] corrected minor typos and employed 
the modern MS renormalization scheme and the ACOT treatment [21] of amplitudes with 
massive quarks (m SjC ^ 0). Very recently, the NLO charm production contributions to 
the full set of electroweak structure functions were calculated [22], including terms that 
are suppressed by m 2 /ME V . In order to apply detector acceptance corrections to the data 
[1], differential NLO distributions were calculated in [23] and [24] that provide the charm 
hadron (D meson) kinematics in terms of the fragmentation variable z and rapidity 77. The 
da / dxdydzdr] code DISCO [24], written by two of the authors of the present article in 
collaboration with D. Mason of NuTeV, exists as an interface to the NuTeV MC event 
generator. Detailed results can be found in the articles listed above. It suffices to say: (i) 
the NLO calculations all agree; and (ii) for the fixed target kinematics under investigation, 
the NLO corrections to the LO results are modest — no bigger than < 20% (see Fig. 1 in 
[17])- 

As mentioned before, the global fits performed in our study are extensions of the full 
NLO CTEQ6 analysis with the addition of constraints due to neutrino dimuon production. 
For the latter process, we have done extensive studies using either the LO formula, Eq. (11), 
or the NLO treatment of [17], Eq. (15). The results obtained in the two cases are quite 
similar. For definiteness, the main results presented in Sec. 4.3 are those obtained by using 
Eq. (11), since the acceptance corrections made to the data set are currently based on a LO 
model, as mentioned earlier. Since we have determined that the NLO corrections to the hard 
cross section are small (compared to experimental errors, for instance), and since we found 
the uncertainty range of the main result (on [S~] ) is much broader than the difference 
between the central values obtained by using Eq. (11) with or without the corrections in 
Eq. (15) (cf. Sec. 5), this approximation does not affect the outcome of our analysis.-" 

1 q — > q for q = s,d. 

s It is certainly desirable to have the inclusive cross sections corrected for acceptance based on NLO models 
(such as [24]), that can be compared to (15) in a full NLO global analysis. This is under active development 
by a theory (CTEQ)-experiment (CCFR-NuTeV) collaboration. 
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4.2 Procedure 

Our analysis is carried out in several stages. First we must find appropriate starting values 
for the fitting parameters. For this purpose, we implement the following steps. 

• We rerun the CTEQ6M global fit with the added CDHSW inclusive neutrino scattering 
data, keeping all other conditions the same. This intermediate fit is extremely close to 
the CTEQ6M one, since the fit to inclusive DIS data is totally dominated by the high 
statistics neutral current experiments on the one hand, and the CCFR charged current 
experiment on the other. 

• We then fix all of the "conventional" parton parameters to their values in this intermediate 
fit, and fit the complete set of data, including the new dimuon data, by varying only the 
parameters associated with the new degrees of freedom in s~ . We obtain results consistent 
with expectations: 

(i) Most of the data sets used in the previous analysis are not affected at all by the 
variation in s~. 

(ii) A few fully inclusive cross sections are slightly affected by the variation of s~ , mainly: 

a. F 3 which depends on u — u + d — d + s — s . . . . 

b. The W ± charge asymmetry which receives a contributions from gs —> W~c. 
These sensitivities to s~ are weak. 

(iii) The CCFR-NuTeV dimuon data sets are the most constraining ones for fitting s"(x). 

We obtain good fits using either the 3-parameter (a, b, c) or the 4- or 5-parameter (a, b, c, d, e) 
versions of Eqs.(8,9). There are not enough constraints to choose among these fits. The 
higher-order polynomials allow oscillatory behavior of s~(x) which the 3-parameter form 
does not. We consider the number of crossings a distinctive property of the physical 
asymmetry s~(x) rather than a volatile function of the continuous fit parameter space. 
As explained in Section 2, one zero-crossing is unavoidably enforced by the sum rule in 
Eq. (5). We are not aware of any solid theoretical argument that would suggest a second 
crossing; nor do we find that the fits show any preference for more than one crossing. We 
therefore restrict the search for best fits in this section to one crossing. This choice also 
seems to be a stable feature of models [2] based on baryon-meson fluctuations. 

• Using these candidate fits as a basis, we perform a second round of fitting allowing the pa- 
rameters associated with s + , Eq.(7), to vary in addition to the s~ variables. This improves 
the fit to all data sets slightly. We observe that the shape of s + (x) now deviates from the 
starting configuration in which s + (x) was set proportional to u(x) + d(x). Defining, as in 
Eq. (1), the strangeness suppression parameter n as the ratio of the momentum fraction 
carried by the strange quarks, [S + ], to that carried by u + d at Qq = 1.3 GeV, we find that 
k may vary in the range 0.3-0.5: \ 2 nas a shallow minimum around k = 0.4. This value 
agrees with previous analyses [11]. 
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Because the experimental constraints are not sufficient to uniquely determine all the s~ 
and s + parameters, we categorize several classes of equally good solutions based on the 
behavior of s~(x)/s + (x) as x — > or x — > 1. 

• We finalize these classes of solutions by allowing all parton parameters to vary so that the 
non-strange parton distributions can adjust themselves to yield the best fit to all the ex- 
perimental data sets. (As one would expect, these final adjustments are generally small.) 
The differences in the x 2 values between the various categories of solutions are not signif- 
icant; i.e. we find nearly degenerate minima with distinctively different s~(x) solutions, 
classified according to their small x behaviour. These solutions do not correspond to iso- 
lated local minima in x 2 space; rather they are to be thought of as specific examples of 
a class of acceptable fits that lie along a nearly flat "valley" along which x 2 changes very 
slowly. 

4.3 Central Results 

The quality of the fits to the global data sets other than the CCFR-NuTeV dimuon data 
remains unaltered from the previous CTEQ6M analysis, so we focus our discussion on the 
strangeness sector. Specifically, we examine closely the asymmetry functions s~(x), S~(x) 
and the momentum integral [S~]. The asymmetry functions from three typical good fits, 
with different behaviors at small x (labeled as classes A,B,C), were previewed in Fig. 2 as 
illustrations. 

In the accompanying table, for each sample fit we list the small-x exponent (3- [s~(x) ~ 
x 13 " , cf. Eq. 10], the integrated momentum fraction [S~], and the relative x 2 values, normal- 
ized to the x 2 °f solution B (xb), which we use as the reference for comparison purposes. 
(Under column "B", we give the absolute x 2 ' s i n parentheses.) 15 To gain some insight on 
the constraints on the strangeness sector due to the various types of experiments, we show 
separately the x 2 values for the dimuon data sets, the inclusive data sets (I) that are ex- 
pected to be somewhat sensitive to s~ (consisting of the CCFR and CDHSW F 3 (x,Q) and 
the CDF H^-lepton asymmetry measurements), and the remaining data sets (II) that are 
only indirectly affected by s~ (the rest of the inclusive data sets). 

k The x 2 values of the dimuon data sets, like those of some other data sets, do not carry rigorous statistical 
significance, because the correlated systematic errors are not available and, hence, cannot be included. In 
the global analysis context, the % 2 value is nevertheless used as the only practical "figure of merit" for the 
fit. The relatively small value of the total x 2 for the dimuon data sets, compared to the number of data 
points, underlines this fact. Under this circumstance, it is common practice to use the normalized y 2 values 
to compare the quality of different fits. 
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#pts 


B+ 


A 


B 


C 


B- 




- 


-0.78 


-0.99 


-0.78 





-0.78 


[S~] x 100 




0.540 


0.312 


0.160 


0.103 


-0.177 


Dimuon 


174 


1.30 


1.02 


1.00 (126) 


1.01 


1.26 


Inclusive I 


194 


0.98 


0.97 


1.00 (141) 


1.03 


1.09 


Inclusive II 


2097 


1.00 


1.00 


1.00 (2349) 


1.00 


1.00 



Table 1. The representative parton distribution sets, arranged in order by the value of [S ]. 

Focusing on the three good fits {A, B, C} first, we note the following features: 

• All three solutions {A,B,C} feature positive [S~\; and the more singular the behavior 
of s~(x) as x — > 0, the higher the value of [S~]. These are natural consequences of 
the strangeness sum rule (equal +/— areas under the curve of s~(x)) and the small x 
suppression of the momentum integral, as discussed earlier in Sec. 2. 

• Solution B is slightly favored over the other two. This, plus the fact that its small- x 
behavior lies in the middle of the favored range, motivates its use as the reference fit. 

• We chose these examples among fits with the simplest parametrizations: all cross the x 
axis only once. With 4- or 5-parameters, which can allow more than one crossing point, 
many solutions can be found that entail oscillatory s~(x). But since the \ 2 values are 
essentially the same as for the simple case, we deem it premature to dwell on complicated 
behaviors, which may be mere artifacts of the parametrization rather than reflections of 
physical constraints. Further studies described in Sec. 4.5 reinforce this point. 

To show how these fits compare with data, we plot in Fig. 3 the ratio of data/theory 
for the reference fit B. The four graphs correspond to the CCFR and NuTeV neutrino and 
antineutrino data sets respectively. The data points are sorted in x-bins, and within each 
x-bin, by y value. We see that the quality of the fit is good, within the experimental 
uncertainties. There are no significant systematic deviations. (The CCFR antineutrino data 
set may appear to be systematically higher than theory. However, upon closer inspection the 
difference is not significant. The data points that lie above theory consist mostly of points 
with large error bars, which tend to catch the attention of the eye; whereas the fit is actually 
dominated by points with small errors, which closely bracket the theory line on both sides. 1 
The value of x 2 /N for this data set is less than 1, comparable to those for the other sets.) 

The parameters for all the fits described in this section are given in detail in the Ap- 
pendix. As is already obvious from Fig. 2, solutions with nearly degenerate \ 2 may corre- 
spond to parametrizations of s~(x) with quite different parameter values, so that a simple 
linear error analysis cannot be applied. This reflects the fact that s~(x) is not well deter- 
mined as a detailed function of x, even when the dimuon data are included in the fit. On 
the other hand, reducing the parameter space to even fewer parameters than our minimal 
set would risk introducing artifacts of an inflexible parametrization. In the next section we 

'This becomes apparent if the data points are re-plotted ordered by the size of the error bars. 
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Data Point Number Data Point Number 

Figure 3: Comparison of data to fit B. The data points are sorted in rr-bins, and within each 
x-bin, by y value. 



will, therefore, apply the Lagrangian multiplier method to deduce the integrated momentum 
asymmetry [S~] and its uncertainty. 



4.4 Range of [S } by the Lagrange Multiplier Method 

Beyond the best fits (A, B, C), we can study the range of [S~] consistent with our global 
analysis in a quantitative way by applying the Lagrange Multiplier (LM) method devel- 
oped in [25]. By varying the Lagrange multiplier parameter, this method explores the entire 
strangeness parameter space in search of solutions with specified values of [S~], i.e., con- 
strained fits. The B~ solution listed in Table 1 was obtained by forcing [S~] = —0.0018 (a 
relatively large negative value, but not as large as the value —0.0027 quoted by [6, 12]). The 
B + solution was generated by forcing [S~] to go in the other (positive) direction until the 
increment of the overall x 2 became comparable to that of B~; this results in [S~] = 0.0054. 

We see from the relevant entries in Table 1 that: (i) the x 2 values of the dimuon data 
sets increase by about 30% in both B ± fits; (ii) the "inclusive I" data sets disfavor the 
negative [S~\; and (iii) the "inclusive II" data sets are completely neutral. These results are 
shown graphically in Fig. 4a, where the square points represent the (relative) x 2 values of 
the dimuon data sets, and the triangle points of the "inclusive I" data sets. (Not shown are 
those for the "inclusive II" data sets, which remain flat (at 1.00).) The LM fits are chosen 
from a large number of fits spanning the entire strangeness parameter space. The pattern of 
dependence of the x 2 values for the dimuon data sets on the value of [S~] is nearly parabolic. 
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This is clear evidence that the dimuon measurement is indeed sensitive to the strangeness 
asymmetry as expected. Further discussion of this observation, including the contrast to the 
sensitivity of other experiments, will be given in Sec. 4.5. 




Figure 4: Values of x 2 /Xb versus [S~] by the Lagrange multiplier method in global analysis. 
Figure (a) show the dimuon data sets (■) and the "Inclusive I" data sets (A) separately, 
and figure (b) shows these data sets combined (■). (Cf. text and Table 1) 

We see from Fig. 4a that, in this series of fits, the dimuon data sets favor a range of 
[S~] centered around 0.0017, whereas the "inclusive I" data sets disfavor negative values of 
[S~]. Fig. 4b shows the dependence of the combined \ 2 °f the two categories of data sets on 
[S~]. (The x 2 °f the remaining data sets used in the global analysis are totally insensitive 
to [S~], cf. Table 1, hence is not included in this plot.) We would like to determine a 
"range of uncertainty" of [S~] from these results. This is far from straightforward because 
of well-known problems shared by all error assessments in global analysis (mainly due to the 
unquantified systematic errors that show up as a lack of statistical compatibility among the 
input data sets). m 

One naive method is to apply the A% 2 = 1 criterion. From the parabola in Fig. 4b, 
which comes from 368 data points, this "estimation-of-parameters criterion" corresponds to 
an uncertainty of [S~] of ±0.0005. (Cf. the lowest horizontal line in Fig. 4b.) It has been 
known, however, that the A% 2 = 1 criterion is unrealistic in global analysis when combining 
data sets with diverse systematic errors from many different experiments [13,25,26]; in this 
circumstance, the overall \ 2 function provides a simple measure of relative goodness-of- 
fit in the minimization process, but it does not have the strict statistical significance of a 
pure parameter fitting problem as presented in textbook examples. This estimate of the 
uncertainty of [S~] is far too small. 

m These difficulties, and practical methods to handle them, are discussed in detail in [13,25,26]. 
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Another often- used method to evaluate "goodness-of-fit" is to apply the cumulative dis- 
tribution function P for the x 2 distribution. One considers unacceptable values of x 2 greater 
than xls ( or X90) where P(x 2 < Xf) = /■ F° r 386 data points, the 68% (90%) criterion 
corresponds to X 2 /Xmin = 1-033 (1.1) respectively. These two criteria are represented by the 
two upper horizontal lines in Fig. 4b. The uncertainty range of [S~] for these two cases are 
±0.002 (0.003) respectively. 

The extensive studies on quantifying uncertainties in the global analysis context [13, 25, 
26] suggest that for this case, a realistic range should be somewhere between the two extreme 
cases shown in Fig. 4b. Hence we adopt the uncertainty range < [S~~] < 0.004 by this 
analysis, which corresponds to the middle horizontal line in the graph (or the xls criterion). 
Whereas a very small strangeness asymmetry, consistent with zero, is not ruled out by this 
criterion, large negative values of [S~] (such as -0.0027, cited in [12]) are strongly disfavored; 
cf. also Table 1. Additional sources of uncertainty will be discussed in the following section. 

4.5 Additional Sources of Uncertainty 

We have performed three series of studies to further assess the rebustness of our main results. 
These will help us to determine a better estimate of the overall uncertainty of [S~]. 

Pure Leading Order Fits Since the experimental analyses of the CCFR-NuTeV dimuon 
data have been done in LO QCD [1, 12], we have carried out a whole series of purely LO 
global analyses, following the same procedures as outline above, in order to provide a basis 
for comparison. The results can be summarized as follows. 

• The overall x 2 f° r the global fit increased by ~ 200 over the comparable fits described 
above; while the x 2 's for the dimuon data sets actually decreased slightly. This is not 
surprising, since the current state of global analysis, with precision data from many ex- 
periments, requires the use of NLO QCD theory - in particular for the collider data with 
typically large perturbative corrections. On the other hand, the new dimuon data still have 
comparably large experimental errors and the NLO corrections are small [0(< 20%)], such 
that an LO fit is adequate for them. 

• We explored the allowed range of strangeness asymmetry [S~] in this LO study under 
different assumptions on the x — > and x — > 1 behavior of the s + (x) and s~(x) functions. 
First, we found that the dependence of Xdimuon 011 * s generally parabolic, rather similar 
to Fig. 4. The width of the distribution is comparable to Fig. 4. The central value for [S~] 
is within the range (0, 0.0015); the exact value depends on the x — > and x — > 1 behavior 
of s ± (x) assumed. 

• We also found that the xLciusive 1 vs - curve, while generally flatter, does "flop around" 
enough for the cases studied so that no clear pattern can be discerned. The specific shape 
of this curve shown in Fig. 3 is not a common characteristic of these fits. 
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Charm Mass Dependence The CCFR-NuTeV dimuon analysis treated the charm mass 
as one of the fit parameters. Their analyses favored a rather high value of m c = 1.6 GeV 
(compared to, e.g., the PDG estimate of 1.0 GeV < m c < 1.4 GeV). The CTEQ global 
analyses are usually done with a fixed value of m c = 1.3 GeV. To see whether the comparison 
between our results is strongly influenced by the choice of the charm mass, we have performed 
several series of fits with m c varying from 1.3 GeV to 1.7 GeV. Again, the general features, 
as described above, stay the same. The central value of [S~] does vary with the choice of 
m c within a given series of fits, but the pattern is not universal. The range over which the 
central value wanders is of the order ~ 0.0015, comparable to the width of the parabola in 
Fig. 4. Unlike the specific analysis of CCFR-NuTeV, the overall \ 2 f° r the global analysis 
does favor a lower value of m c . 

Dependence on Decay and Fragmentation Model To estimate the dependence of 
our results on the model used to convert the measured dimuon cross sections to structure 
functions for charm production, we repeated our analyses using an alternative conversion 
table provided by the CCFR-NuTeV collaboration. 11 This alternative table is based on 
Buras-Gaemers PDFs used in CCFR-NuTeV analyses with Collins-Spiller fragmentation 
functions. It is similarly tuned to detailed features of the measured dimuon cross sections as 
that described in Sec. 4.2.° The results obtained from the alternative fits are, again, similar 
to those described earlier. The Xdimuon vs - [&~] parabola generally has the same width as in 
Fig. 4. The central value for [S~] is in the range < [S~~] < 0.0015 — in the lower half of the 
range quoted at the end of Sec. 4.4. The dependence of xfnciusive i on shows no definitive 
trend. 

Taken together, the results of the additional studies described in these three paragraphs lead 
to several conclusions, (i) The general features described in Sees. 4.3 and 4.4 are robust; (ii) 
The central value of [S~] wanders around within a range that is consistent with the width of 
the Xdimuon vs - parabola; and (iii) These additional results do not significantly change 
the estimates of the previous section, except to shift the central estimated value of [S~] to 
a slightly lower value, and to extend the range of uncertainty on the lower side somewhat. 
The envelope of these additional uncertainties provides an estimated range of uncertainty of 
the strangeness asymmetry of 9 

-0.001 < [S-] < 0.004 . (16) 

This large range reflects both the limit of current experimental constraints and the consider- 
able theoretical uncertainty, as explicitly discussed in the text. The theoretical uncertainties 
"We thank Kevin MacFarland for supplying this table. 

°However, since our CTEQ6-likc PDFs arc rather different from the CCFR Buras-Gaemers PDFs, it is 
not clear how good the approximation is to use this conversion table. That is, the self-consistency of the 
procedure is not assured. 

p The uncertainties from the various sources need not be combined in quadrature, because they are not 
statistically independent sources, but rather systematic uncertainties of the theory. 
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can be reduced in a refined NLO analysis; the results remain to be seen. The limitations on 
the experimental constraints will remain, until new experiments are done . 

5 Comparisons to previous studies 

A comprehensive global QCD analysis with emphasis on the strangeness sector has been 
carried out previously by BPZ [4]. q Without the dimuon data, which are directly sensitive to 
strangeness, the results of BPZ implicitly rely on small differences between large neutral- and 
charged- current inclusive cross sections from different experiments. The latest representative 
s~(x) and S~(x) functions extracted by BPZ are shown in Fig. 1, along with our reference fit 
B. The main feature of the BPZ curves is a positive bump at rather large x. r This feature has 
been attributed to the influence of the CDHSW data, particularly when re-analyzed at the 
cross section level along with the other DIS experiments. Their conclusion that data favor 
a positive value of the momentum integral [S~] is in general agreement with our detailed 
study based on the LM method. However, the different shapes of s~(x) seen in Fig. 1 clearly 
underline the difference in inputs: (i) our results are mainly dictated by the CCFR-NuTeV 
dimuon data (which are not present in the BPZ analysis); (ii) their results rely on a delicate 
analysis of DIS cross-section data (not matched in our structure function analysis); and (iii) 
the difference in flexibility of the parametrizations of the non-perturbative input functions 
can influence the results. 

The CCFR and NuTeV collaborations performed separate and combined analyses of s 
and s [1], based on their own dimuon and inclusive cross sections. To parameterize the s(x) 
and s(x) distributions, they chose the model formula 

/ s(x, Q) \ _ u(x, Q) + djx, Q) ( K (l - x) a \ 

V K*,Q) J 2 V «(l-x) a J [ } 

for all {x, Q}, where k, k, a, a are fitting parameters. 

Curves representing the general behavior of the model (17) at Q 2 = lOGeV 2 , with 
parameter values (k, k, a, a) taken from [1], have been shown in Fig. 1 for comparison with 
the other distributions. While this model might be acceptable for a limited range of x and 
Q, it leads to problems when extrapolated to general {x, Q} values: (i) the strangeness 
number sum rule [s~] = is violated (in fact, the integral [s~] diverges unless 1% — K^j j clS 
is the momentum sum rule, Eq. 6; (ii) the QCD evolution equations are violated by the 

q As mentioned in Sec. 4, BPZ work directly with DIS cross sections (instead of structure functions), with 
detailed attention to systematic errors and other sources of uncertainties. 

r The BPZ curve displayed includes the more recent analysis "with CCFR (inclusive data)." The original 
"without CCFR" solution in [4] has a more pronounced large-x bump and a smaller negative region. An 
s~ (x) function of such magnitude at large x is in disagreement with the CCFR-NuTeV analysis; in addition, 
it would make s{x) and s(x) behave quite differently compared with the non-strange sea quarks and the 
gluon, i.e., ~ (1 — x) p , with p in the range 5 — 10. 
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Q-dependent paramaterizations of the PDFs. s The first problem is evident in Fig. 1.* 

It could be argued that, since the experiment only covers a limited range of x, the 
enforcement of the sum rules is not critical in extracting limited information on s(x, Q) 
and s(x,Q). Given this, the CCFR curve in Fig. 1 implies negative s~(x) over most of the 
experimental x range (0.01 - 0.3). u The uncertainty is smaller at the lower x end because 
of better statistics. It is this feature of the data that we invoked in the general discussion of 
Sec. 2. 

To describe the behavior of s~(x) over the full x range, the strangeness number sum 
rule must be enforced. The sum rule also provides a powerful theoretical constraint on the 
data analysis, as demonstrated in Sees. 2 and 4. For functions to be candidate universal 
parton distributions of the PQCD formalism, they must satisfy both the sum rules and the 
QCD evolution equations. 

6 Conclusion 

We find a range of solutions in the strangeness sector that are consistent with all rele- 
vant world data used in the global analysis. The dimuon data are vital in constraining the 
strangeness asymmetry parameters. The constraints provided by other inclusive measure- 
ments, labeled as "inclusive I" in the text, are consistent with those provided by dimuon 
data, although much weaker. The allowed solutions generally prefer the momentum integral 
[S~] = Jq x[s(x) — s(x)) dx to be positive. This conclusion is quite robust, and it follows 
from the basic properties of PQCD and from qualitative features of the experimental data. 
However, the size of this strangeness momentum asymmetry is still quite uncertain; we can 
only estimate that [S~] lies in the range from —0.001 to +0.004. The Lagrange Multiplier 
method explicitly demonstrates that both the dimuon data and the "inclusive I" data sets 
strongly disfavor a large negative value of [S~], although they may still be consistent with 
zero asymmetry. 

The fact that [S~] has a large uncertainty has significant implications for the precision 
measurement of the weak mixing angle, sin 2 9w, from neutrino scattering. This issue is 
studied separately in Ref. [8]. 

This paper marks the first global QCD analysis incorporating direct experimental con- 
straints on the strangeness sector. We have so far focused only on the strangeness asymmetry, 
which represents a new frontier in parton degrees of freedom. Much still needs to be done 
to improve the treatment of the dimuon data (to true NLO accuracy), and to fully explore 

s These problems are independent of the specific issues of the strange quark asymmetry. 

'Re-analysis of the CCFR-NuTeV data by the experimental group, taking into account these issues, are 
underway Initial results from partial implementations of the above-mentioned theoretical constraints were 
reported by P. Spentzouris at International Workshop on Weak Interactions and Neutrinos 2003 (WIN03), 
Lake Geneva, October, 2003. 

"This is consistent with the fact that (without the constraint of sum rules) [12] quotes [5—] = — 0.0027 ± 
0.0013. 
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all the degrees of freedom associated with s + (x) and s~(x). As progress is made on these 
fronts, the uncertainty on [S~] will no doubt decrease as well. 

Note added: After this manuscript was completed an investigation of the 3-loop perturba- 
tive strangeness asymmetry was presented in [27]. 

Acknowledgment 

We thank members of the CCFR and NuTeV collaboration, in particular T. Bolton and 
M. Goncharov, for discussions about the dimuon data, and assistance in their use, and Kevin 
MacFarland, Dave Mason, and Panagiotis Spentzouris for useful comments and suggestions. 
We also thank Benjamin Portheault for discussions of the BPZ work and for assistance in 
generating their results for comparison. F.O. acknowledges the hospitality of MSU and BNL 
where a portion of this work was performed. This research was supported by the National 
Science Foundation (grant No. 0100677), RIKEN, Brookhaven National Laboratory, and 
the U.S. Department of Energy (Contract No. DE-AC02-98CH10886, and No. DE-FG03- 
95ER40908), and by the Lightner-Sams Foundation. 

Appendix 

In this appendix, we provide some detailed information on the parametrization of the non- 
perturbative parton distribution functions used in our global analysis work, as well as the 
values of the parameters for the representative fits. 

The non-strange parton distributions at Q — Qo = 1.3 GeV are parametrized as in 
CTEQ6. We use 

x f( Xi Q ) = A x Al (1 - x) M e A3X (1 + e M x) As (18) 

with independent parameters for parton flavor combinations u v = u — u, d v = d — d, g, and 
u + d . To distinguish the d and u distributions, we parametrize the ratio d/u, as a sum of 
two terms: 

d(x, Q )/u(x, Qo) = At x M (1 - x) As + (1 + A 4 x) (1 - x) M . (19) 

The strangeness sector is parametrized according to the description of Sec. (3) with 

xs + (x, Qo) = A x Al (1 - x) M e Mx (1 + e M x) M 

where the 2,3,4,5 coefficients are either equated to those of u + d, or allowed to vary 
independently, depending on the particular fit being performed; and 

s~(x,Qq) = s + (x,Q )tanh[A x Al (l-x) M (l - -^J (l + A A x + A 5 x 2 )} 

As an example, the "standard fit" (B-fit) described in Sec. 4. 3 has the following coeffi- 
cients. 
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B-fit 


A 


A x 


A 2 


^3 


A, 


A 5 


d v 


1.55891 


0.62792 


5.08865 


0.70688 


-0.26338 


3.00000 


u v 


1.73352 


0.55260 


2.90090 


-2.63846 


1.45774 


1.85987 


9 


31.24912 


0.52014 


2.38230 


4.23010 


2.33765 


-3.00000 


d/u 




10.23811 


5.19599 


14.85860 


17.00000 


8.68941 


u + d 


0.06758 


-0.29681 


7.71700 


-0.82223 


4.43481 


0.66711 


s + s 


0.04966 


0.03510 


7.44149 


-1.44570 


5.13400 


0.59659 


s — s 


-0.08438 


0.18803 


2.36708 


0.04590 


0.00000 


0.00000 



And the other sample fits have coefficients: 



A-fit 


A 


Ax 


A 2 


^3 


A, 


A 5 


d v 


1.55485 


0.62740 


5.09550 


0.70010 


-0.25130 


3.00000 


u v 


1.73345 


0.55260 


2.90090 


-2.63740 


1.45810 


1.85910 


9 


31.05518 


0.51840 


2.38230 


4.24730 


-2.33790 


3.00000 


d/u 




10.23920 


5.19690 


14.85860 


17.00000 


8.68530 


u + d 


0.06677 


-0.29810 


7.71700 


-0.80540 


4.45630 


0.66420 


s + s 


0.05133 


0.06567 


7.59880 


-1.43972 


5.13400 


0.61715 


s — s 


-0.02000 


-0.06562 


-0.48403 


0.01987 


-0.58837 


0.00000 



C-fit 


A 


Ax 


A 2 


^3 


A, 


A 5 


d v 


1.39587 


0.60594 


4.75379 


-0.75291 


0.31387 


3.00000 


u v 


1.72080 


0.55260 


2.90090 


-2.36331 


1.60950 


1.59812 


9 


29.40547 


0.50423 


2.38230 


4.41409 


2.34266 


-3.00000 


d/u 




10.23202 


5.19618 


14.85860 


17.00000 


8.62620 


u + d 


0.06479 


-0.30038 


7.71700 


-0.68604 


4.61325 


0.62988 


s + s 


0.04289 


0.00809 


7.71700 


-1.22365 


5.13400 


0.62988 


s — s 


-2.36926 


0.99191 


8.23499 


0.07884 


0.00000 


0.00000 



B+-fit 


A 


Ax 


A 2 


^3 


A, 


A 5 


d v 


1.43611 


0.61170 


4.73270 


-0.67420 


0.24920 


3.00000 


u v 


1.71921 


0.55260 


2.90090 


-2.39220 


1.60430 


1.61490 


9 


29.76781 


0.50800 


2.38230 


4.35570 


-2.33860 


3.00000 


d/u 




10.19560 


5.16810 


14.85860 


17.00000 


8.69760 


u + d 


0.06729 


-0.29650 


7.71700 


-0.75670 


4.52290 


0.64380 


s + s 


0.03456 


0.00210 


8.23420 


-1.26970 


5.13410 


0.72501 


s — s 


-0.40480 


0.22103 


3.40190 


0.04701 


0.31550 


0.00000 



20 



B_-fit 


A 


Ai 


A 2 


^3 


A 4 


A 5 


d v 


1.43611 


0.61170 


4.73270 


-0.67420 


0.24920 


3.00000 


u v 


1.71921 


0.55260 


2.90090 


-2.39220 


1.60430 


1.61490 


9 


29.76781 


0.50800 


2.38230 


4.35570 


-2.33860 


3.00000 


d/u 




10.19560 


5.16810 


14.85860 


17.00000 


8.69760 


u + d 


0.06717 


-0.29650 


7.71700 


-0.75670 


4.52290 


0.64380 


s + s 


0.04356 


0.00210 


7.33918 


-1.26970 


5.13400 


0.60083 


s — s 


0.01781 


0.22103 


-15.02691 


0.22693 


-1.23666 


0.00000 
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